GeoMx_EDA
Jiwoon Park
I4 GeoMx Skin spatial transcriptomics analysis
The analysis in this markdown starts from the Q3-normalized count matrices and metadata, geomx (RDS file available via NASA GeneLab repository: URL to be updated). We also provide estimated cell proportion object geomx_cellprop, deconvoluted from the HCA reference dataset (source link to be updated).
source("../Code/HelperFunctions.R")## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.geomx <- readRDS("../Data/GeoMx_merged_pub.RDS")
geomx_cellprop <- readRDS("../Data/i4GeoMxData_cellproportionestimates.RDS")From initial clustering analysis, two outlier samples were detected from the PCA analysis: (1) C001_Pre_OE_1 and (2) C004_Post_VA_1. the final object already excludes those two samples. To generate the UMAP plot in the manuscript:
# PCA
pca <- prcomp((geomx[,35:ncol(geomx)]), scale=TRUE)
pcares <- summary(pca)
d <- data.frame(x = pca$x[,1], y = pca$x[,2], sample = geomx$ROInameJP, Timepoint = geomx$Timepoint, tissue.types = geomx$ROIType, crew = geomx$SlideName)
ggplot(d, aes(x, y, shape = Timepoint, color = tissue.types)) +
geom_point(size = 3) +
theme_pub +
xlab(paste("PC1: ", round(pcares$importance[2,1],3)*100, "% variance", sep="")) +
ylab(paste("PC2: ", round(pcares$importance[2,2],3)*100, "% variance", sep=""))
# UMAP
umap.rld = umap::umap(geomx[,35:ncol(geomx)])
umap.df = umap.rld$layout %>% as.data.frame()
umap.df$sample = geomx$ROInameJP; umap.df$Timepoint = geomx$Timepoint; umap.df$Type = geomx$ROIType; umap.df$crew = geomx$SlideName
ggplot(umap.df, aes(V1, V2, shape = Timepoint, color = Type)) +
geom_point(size = 4, alpha = 0.7) +
theme_bw() +
xlab("UMAP1") +
ylab("UMAP2") +
stat_ellipse(linetype=2, aes(group=Type)) +
facet_wrap(~Timepoint) +
scale_color_brewer(palette = "Dark2")## Too few points to calculate an ellipse## Warning: Removed 1 row containing missing values (`geom_path()`).
Differential expression analysis
# (1) Pre- and Post- spaceflight comparisons
ddsMatrix <- DESeq2::DESeqDataSetFromMatrix(countData=geomx[,35:ncol(geomx)] %>% t() %>% round(), colData=geomx[,1:35], design=~Timepoint)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factorsdds <- DESeq(ddsMatrix)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 102 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingvsd <- vst(dds)## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.overall <- results(dds, contrast=c("Timepoint", "Post", "Pre")) %>% as.data.frame()
overall_sig = overall %>% filter(padj < 0.05 & log2FoldChange > 0)
# (2) Pre- and Post- spaceflight comparisons, further divided by tissue regions
ddsMatrix <- DESeq2::DESeqDataSetFromMatrix(countData=geomx[,35:ncol(geomx)] %>% t() %>% round(), colData=geomx[,1:35], design=~ROITime)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factorsdds <- DESeq(ddsMatrix)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 37 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDEG_OE <- results(dds, contrast=c("ROITime", "Post_OE", "Pre_OE")) %>% as.data.frame()
DEG_IE <- results(dds, contrast=c("ROITime", "Post_IE", "Pre_IE")) %>% as.data.frame()
DEG_OD <- results(dds, contrast=c("ROITime", "Post_OD", "Pre_OD")) %>% as.data.frame()
DEG_VA <- results(dds, contrast=c("ROITime", "Post_VA", "Pre_VA")) %>% as.data.frame()
padj_cutoff=0.5; log2FoldChange_cutoff=0.1
DEG_OE_sig <- DEG_OE %>% filter(padj < padj_cutoff & abs(log2FoldChange) > log2FoldChange_cutoff)
DEG_IE_sig <- DEG_IE %>% filter(padj < padj_cutoff & abs(log2FoldChange) > log2FoldChange_cutoff)
DEG_OD_sig <- DEG_OD %>% filter(padj < padj_cutoff & abs(log2FoldChange) > log2FoldChange_cutoff)
DEG_VA_sig <- DEG_VA %>% filter(padj < padj_cutoff & abs(log2FoldChange) > log2FoldChange_cutoff)
DEG_VA_sig <- DEG_VA %>% filter(pvalue < 0.01 & abs(log2FoldChange) > 0)
# compare DEGs
venndf <- list(overall = rownames(overall_sig),
OE = rownames(DEG_OE_sig),
IE = rownames(DEG_IE_sig),
OD = rownames(DEG_OD_sig),
VA = rownames(DEG_VA_sig))
#plot(euler(venndf, shape="ellipse"), quantities=TRUE)
UpSetR::upset(fromList(venndf), keep.order=TRUE, order.by = "freq")
FC_cutoff = 0.5; adjp_cutoff = 0.1
# get pathways to use
pathways.hallmark <- gmtPathways("../../Codes/pathwaydata/h.all.v7.5.1.symbols.gmt")
pathways.c2 <- gmtPathways("../../Codes/pathwaydata/c2.all.v7.5.1.symbols.gmt")
pathways.c5 <- gmtPathways("../../Codes/pathwaydata/c5.all.v7.5.1.symbols.gmt")
pathways.c7 <- gmtPathways("../../Codes/pathwaydata/c7.all.v7.5.1.symbols.gmt")
pathwaylist <- list("hallmark" = pathways.hallmark, "canonical" = pathways.c2, "gene_ontology" = pathways.c5, "immune" = pathways.c7)
# generate volcano plots: overall comparison
genedf <- overall
processed_df <- genedf %>% process_voldf()
subdf <- processed_df[[2]] %>% as.data.frame()
subdf_degs <- c(nrow(subdf[subdf$direction == "up",]), nrow(subdf[subdf$direction == "down",]))
volcano_plot(df = processed_df[[1]], subdf = processed_df[[2]],
plotname = paste("Volcano plot: ", "overall Post vs. Pre flight", ", |FC| >", FC_cutoff, " and adj.p-value <", adjp_cutoff, sep="")) +
labs(subtitle = paste("(", subdf_degs[1], " up-regulated and ", subdf_degs[2], " down-regulated genes.)", sep="")) +
theme(plot.subtitle = element_text(hjust = 1))## Warning: ggrepel: 41 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# or group of volcano plots by region types
ggarrange(run_volcano(DEG_OE, "OE"), run_volcano(DEG_IE, "IE"), run_volcano(DEG_OD, "OD"), run_volcano(DEG_VA, "VA"), ncol=2, nrow=2)## Warning: Removed 496 rows containing missing values (`geom_label_repel()`).## Warning: Removed 350 rows containing missing values (`geom_label_repel()`).## Warning: ggrepel: 109 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 105 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 445 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# run pathway analysis
# can be looped to run all analyses
for(path in c("hallmark", "canonical", "gene_ontology", "immune")){
pathway_db <- path
pathway_choice <- pathwaylist[[pathway_db]]
pathway_df <- genedf
pathway_df$GENENAME <- rownames(pathway_df)
#can be saved with: write.csv(pathway_analysis(df = pathway_df , pathway = pathway_choice), paste("../Analysis/PathwayAnalysis/Pathway_All_", path, "_VA_all_DEGs.csv", sep=""))
}
# or run individually
genedf <- overall
pathway_db <- "hallmark"
pathway_choice <- pathwaylist[[pathway_db]]
pathway_df <- genedf
pathway_df$GENENAME <- rownames(pathway_df)
plot_pathways(fgseaRes = pathway_analysis(df = pathway_df , pathway = pathway_choice), p_cut = 0.05)
# visualize into heatmap or barplot sets for multi-condition comparisons
pathway_db <- "hallmark"
pathway_choice <- pathwaylist[[pathway_db]]
genedf_OE = DEG_OE; genedf_OE$GENENAME <- rownames(genedf_OE); pres_OE = pathway_analysis(df = genedf_OE, pathway = pathway_choice)
genedf_IE = DEG_IE; genedf_IE$GENENAME <- rownames(genedf_IE); pres_IE = pathway_analysis(df = genedf_IE, pathway = pathway_choice)
genedf_OD = DEG_OD; genedf_OD$GENENAME <- rownames(genedf_OD); pres_OD = pathway_analysis(df = genedf_OD, pathway = pathway_choice)
genedf_VA = DEG_VA; genedf_VA$GENENAME <- rownames(genedf_VA); pres_VA = pathway_analysis(df = genedf_VA, pathway = pathway_choice)
pres_OE$ROIType <- "OE"; pres_IE$ROIType <- "IE"; pres_OD$ROIType <- "OD"; pres_VA$ROIType <- "VA"
pres_all = rbind(pres_OE, pres_IE, pres_OD, pres_VA)
pres_all$ROIType <- factor(pres_all$ROIType, levels=c("OE", "IE", "OD", "VA"))
pres_all$pathway = gsub("HALLMARK_", "", pres_all$pathway)
pres_all$pathway = gsub("_", " ", pres_all$pathway)
plot_pathways(pres_all, p_cut = 1) + facet_wrap(~ROIType, ncol=4)
path_heatdf = pres_all[,c("pathway", "NES", "ROIType")] %>% tidyr::spread(ROIType, NES)
rownames(path_heatdf) = path_heatdf$pathway
pheatmap::pheatmap(path_heatdf[,2:5] %>% t(), scale="none", cluster_rows = FALSE)
Below shows multiple methods we used to visualize cellular proportion changes:
ggplot(geomx_cellprop %>% tidyr::gather(celltype, proportions, 35:ncol(geomx_cellprop)) %>% group_by(ROInameJP, celltype),
aes_string(y="proportions", x="Timepoint", fill="Timepoint")) +
geom_boxplot(color="grey30", alpha=0.75) +
#scale_fill_brewer(palette="Set1") +
scale_fill_manual(values=c("#377eb8", "#e41a1c")) +
facet_grid(~celltype) +
ggtitle(paste("Cell Proportion Change during flight, overall", sep=" ")) + stat_compare_means(label="p.signif", hjust=-1) + theme_bw()## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation ideoms with `aes()`
cellprop_plot <- geomx_cellprop %>% tidyr::gather(celltype, proportions, 35:ncol(geomx_cellprop)) %>% group_by(ROInameJP, celltype)
ggplot(cellprop_plot, aes(x=ROInameJP, y=proportions, fill=celltype)) +
geom_bar(stat="identity", position="fill", width=1, color="white") +
#coord_polar("y", start=0) +
#facet_grid(~Timepoint) +
theme_void()
cellprop_plot <- geomx_cellprop %>% tidyr::gather(celltype, proportions, 35:ncol(geomx_cellprop)) %>% group_by(Timepoint, celltype) %>% summarise(AverageProp = mean(proportions))## `summarise()` has grouped output by 'Timepoint'. You can override using the
## `.groups` argument.ggplot(cellprop_plot, aes(x="", y=AverageProp, fill=celltype)) +
geom_bar(stat="identity", position="fill", width=1, color="white") +
#coord_polar("y", start=0) +
facet_wrap(~Timepoint) +
theme_void()
cellprop_plot <- geomx_cellprop %>% tidyr::gather(celltype, proportions, 35:ncol(geomx_cellprop)) %>% group_by(ROIType, Timepoint, celltype) %>% summarise(AverageProp = mean(proportions))## `summarise()` has grouped output by 'ROIType', 'Timepoint'. You can override
## using the `.groups` argument.cellprop_plot$Type = paste(cellprop_plot$ROIType, cellprop_plot$Timepoint, sep="_")
ggplot(cellprop_plot, aes(x="", y=AverageProp, fill=celltype)) +
geom_bar(stat="identity", position="fill", width=1, color="white") +
coord_polar("y", start=0) +
facet_wrap(~Type, ncol=4) + theme_void()
cellprop_plot <- geomx_cellprop %>% tidyr::gather(celltype, proportions, 35:ncol(geomx_cellprop)) %>% group_by(Timepoint, celltype, ROIType) %>% summarise(AverageProp = mean(proportions))## `summarise()` has grouped output by 'Timepoint', 'celltype'. You can override
## using the `.groups` argument.cellprop_plot0 <- cellprop_plot[which(cellprop_plot$Timepoint == "Pre"),]
cellprop_plot$AveragePropFC <- (cellprop_plot$AverageProp+0.000001)/(cellprop_plot0$AverageProp+0.000001)
ggplot(cellprop_plot[which(cellprop_plot$Timepoint == "Post"),], aes(x=Timepoint, y=AveragePropFC, fill=celltype)) +
geom_bar(stat="identity", position="dodge", width=1, color="white") +
#coord_polar("y", start=0) +
ylim(0,2) +
facet_grid(~ROIType) +
theme_bw()## Warning: Removed 1 rows containing missing values (`geom_bar()`).
cellprop_hdf = cellprop_plot[which(cellprop_plot$Timepoint == "Post"),c(2,3,5)] %>% group_by(ROIType) %>% tidyr::spread("celltype", "AveragePropFC") %>% as.data.frame()
rownames(cellprop_hdf) = cellprop_hdf$ROIType
cellprop_hdf = cellprop_hdf[,2:ncol(cellprop_hdf)]
cellprop_hdf = cellprop_hdf[c("OE", "IE", "OD", "VA"),]
pheatmap(as.matrix(cellprop_hdf) %>% scale(), scale="row", cluster_rows = FALSE)
cellprop_hdf[,2:ncol(cellprop_hdf)] %>% as.matrix()## Fibroblast Keratinocyte Lymphatic_EC Macrophages_DC Melanocyte Pericyte
## OE 1.3784230 0.7685123 1.0836908 0.9822108 0.8359297 1.0841828
## IE 0.4785430 0.6638136 0.6768361 0.6561741 0.6309900 0.6071922
## OD 0.5913574 1.0000000 0.7277393 0.7303084 0.6766703 0.6229705
## VA 0.4871110 1.0000000 0.6452769 1.5489680 0.9365917 0.4891048
## T Vascular_EC
## OE 1.1773709 1.000000e+00
## IE 0.5067802 1.000000e+00
## OD 0.4758284 9.222414e-01
## VA 0.2541789 2.964439e+04Correlation plots to see cellular microenvironment changes
b <- as.character(unique(geomx_cellprop$ROITime))
df_forcorplot <- data.frame(geomx_cellprop[,35:43]) #data.frame(geomx_cellprop[,c(19,17,18,15,13,16,11,20,14,8,9,7,12,6)]) # cell type order
corplot_byROItime <- vector('list', length(b))
for(i in 1:length(b)) {
corplot_byROItime[[i]] <- ggcorrplot(cor((df_forcorplot[geomx_cellprop$ROITime == b[i],])),
hc.order = FALSE,
colors = c("blue", "grey", "red"),
insig = "pch", #"blank",
p.mat = cor_pmat((df_forcorplot[geomx_cellprop$ROITime == b[i],]))) +
theme_pubclean() +
ggtitle(paste(b[i])) +
theme(axis.text.x = element_text(angle=45, hjust=1),
axis.title.x = element_blank(),
axis.title.y = element_blank())
}## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## Warning in cor((df_forcorplot[geomx_cellprop$ROITime == b[i], ])): the standard
## deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zeroggarrange(plotlist = corplot_byROItime, nrow = 2, ncol = 4)
Individual expression boxplots can be drawn as follows:
geneplot_df <- geomx[geomx$ROIType %in% c("OE", "IE", "OD", "VA"),]
geneplot_df$ROIType <- factor(geneplot_df$ROIType, levels=c("OE", "IE", "OD", "VA"))
genenames = "DPP4"
ggplot(geneplot_df, aes_string(y=genenames, x="Timepoint", fill="Timepoint")) +
geom_boxplot(color="grey30", alpha=0.75) +
scale_fill_brewer(palette="Set1") +
coord_flip() +
ggtitle(paste("Expression Change during flight:", genenames, sep=" ")) +
facet_wrap(~ROIType, ncol=1) +
stat_compare_means(label="p.signif", hjust=-.7) +
theme_pubr() + theme(legend.position="none")
Expression profiles of multiple genes can be shown in a form of heatmaps:
geneplot_df <- geomx[geomx$ROIType %in% c("OE", "IE", "OD", "VA"),]
geneplot_df$ROIType <- factor(geneplot_df$ROIType, levels=c("OE", "IE", "OD", "VA"))
annoCol <- list(SlideName = c(C001 = "#f0e7e7", C002 = "#c1440e", C003 = "#e77d11", C004 = "#fda600"),
ROIType = c(OE = "#1b9e77", IE = "#d95f02", OD = "#7570b3", VA = "#e7298a"),
Timepoint = c(Pre = "#377eb8", Post = "#e41a1c"))
avgexp = geneplot_df %>% tidyr::gather(genename, value, 35:ncol(geneplot_df)) %>% group_by(ROIType, Timepoint, genename) %>% summarise(AverageExp = mean(value)) %>% tidyr::spread(genename, AverageExp)## `summarise()` has grouped output by 'ROIType', 'Timepoint'. You can override
## using the `.groups` argument.sdexp = geneplot_df %>% tidyr::gather(genename, value, 35:ncol(geneplot_df)) %>% group_by(ROIType, Timepoint, genename) %>% summarise(AverageExp = sd(value)) %>% tidyr::spread(genename, AverageExp)## `summarise()` has grouped output by 'ROIType', 'Timepoint'. You can override
## using the `.groups` argument.rownames(avgexp) = paste(avgexp$Timepoint, avgexp$ROIType, sep="_")## Warning: Setting row names on a tibble is deprecated.avgexp_FC <- avgexp[c(1,3,5,7),3:ncol(avgexp)]/avgexp[c(2,4,6,8),3:ncol(avgexp)]
rownames(avgexp_FC) <- avgexp$ROIType[c(1,3,5,7)]
genenames = c("STAT3", "STAT5B", "JAK1", "FLG", "SPINK5", "DSG1", "DSP", "CDSN", "ADAM17", "EGFR") # can check whether they exist in geomx by the following line: genenames <- genenames[which(genenames %in% colnames(geneplot_df))]
# to plot fold changes
geneplot_ref = colMeans(geneplot_df[geneplot_df$Timepoint=="Pre",35:ncol(geneplot_df)])
pheatmap::pheatmap(geneplot_df[order(geneplot_df$ROIType, geneplot_df$Timepoint, geneplot_df$SlideName),genenames[1:length(genenames)]],# %>% scale(),
color = colorRampPalette(rev(brewer.pal(n = 10, name = "RdGy")))(100),
border_color = "grey60",
cluster_cols = TRUE, cluster_rows = FALSE,
legend = TRUE, scale = "column",
annotation_row = geneplot_df[,c("Timepoint", "ROIType", "SlideName")],
annotation_colors = annoCol,
show_colnames = TRUE, show_rownames = FALSE)
geneplot_df_heatmap <- geneplot_df[geneplot_df$Timepoint=="Post", genenames]/geneplot_ref[genenames]
pheatmap::pheatmap(geneplot_df_heatmap[order(geneplot_df$ROIType[geneplot_df$Timepoint=="Post"], geneplot_df$Timepoint[geneplot_df$Timepoint=="Post"], geneplot_df$SlideName[geneplot_df$Timepoint=="Post"]),],# %>% scale(),
color = colorRampPalette(rev(brewer.pal(n = 10, name = "RdGy")))(100),
border_color = "grey60",
cluster_cols = T, cluster_rows = F,
legend = TRUE, scale = "column",
annotation_row = geneplot_df[,c("ROIType", "SlideName")],
annotation_colors = annoCol,
show_colnames = T, show_rownames = F)
annoCol2 <- list(ROIType = c(OE = "#1b9e77", IE = "#d95f02", OD = "#7570b3", VA = "#e7298a"),
Timepoint = c(Post = "#e41a1c", Pre = "#377eb8"))
annoCol_df <- avgexp[,c("ROIType", "Timepoint")]
annoCol_df$Timepoint <- factor(annoCol_df$Timepoint, levels = c("Pre", "Post"))
avgexp_plot = avgexp[order(avgexp$ROIType, avgexp$Timepoint),c("ROIType", "Timepoint", genenames)]
rownames(avgexp_plot) = paste(avgexp_plot$Timepoint, avgexp_plot$ROIType, sep="_")## Warning: Setting row names on a tibble is deprecated.avgexp_plot0 = avgexp_plot[,3:ncol(avgexp_plot)]; rownames(avgexp_plot0) = rownames(avgexp_plot)## Warning: Setting row names on a tibble is deprecated.pheatmap::pheatmap(avgexp_plot0, #%>% scale(),
color = colorRampPalette(rev(brewer.pal(n = 10, name = "RdGy")))(100),
border_color = "grey60",
cluster_cols = T, cluster_rows = F,
legend = TRUE, scale = "column",
show_colnames = T, show_rownames = T)
avgexp_plot1 <- log2(avgexp_plot0[c(1,3,5,7),])-log2(avgexp_plot0[c(2,4,6,8),]); rownames(avgexp_plot1) <- rownames(avgexp_plot0)[c(1,3,5,7)]
pheatmap::pheatmap(avgexp_plot1, #%>% scale(),
color = colorRampPalette(rev(brewer.pal(n = 10, name = "RdGy")))(100),
border_color = "grey60",
cluster_cols = T, cluster_rows = F,
legend = TRUE, scale = "column",
show_colnames = T, show_rownames = T)
Single sample GSEA (ssGSEA) was performed to see gene enrichments by regions of interest (ROI):
# make genesets
geneset <- readRDS("../Data/Custom_GeneSets.RDS")
ssgsea_rawmat <- geomx[,35:ncol(geomx)] %>% t()
colnames(ssgsea_rawmat) <- geomx$ROInameJP
ssgsea_res <- gsva(ssgsea_rawmat, method="ssgsea", geneset)## Estimating ssGSEA scores for 31 gene sets.
## [1] "Calculating ranks..."
## [1] "Calculating absolute values from ranks..."
##
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|===== | 6%
|
|===== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 35%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 49%
|
|=================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
##
## [1] "Normalizing..."ssgsea_plotdf <- t(ssgsea_res) %>% as.data.frame(); ssgsea_plotdf$ROInameJP <- rownames(ssgsea_plotdf)
ssgsea_plotdf01 <- dplyr::inner_join(geomx[,1:35], ssgsea_plotdf, by="ROInameJP")
ssgsea_plotdf01 <- ssgsea_plotdf01[ssgsea_plotdf01$ROIType %in% c("OE", "IE", "OD", "VA"),]
ssgsea_plotdf01$ROIType <- factor(ssgsea_plotdf01$ROIType, levels=c("OE", "IE", "OD", "VA"))
pheatmap((ssgsea_plotdf01[order(ssgsea_plotdf01$ROIType, ssgsea_plotdf01$Timepoint, ssgsea_plotdf01$SlideName),36:ncol(ssgsea_plotdf01)]) %>% t(),
annotation_col = ssgsea_plotdf01[,c("Timepoint", "ROIType", "SlideName")],
color = colorRampPalette(rev(brewer.pal(n = 10, name = "BrBG")))(100),
annotation_colors = annoCol,
show_rownames = TRUE, show_colnames = FALSE,
scale = "row", cluster_rows = TRUE, cluster_cols = FALSE)
ssgsea_avg = ssgsea_plotdf01 %>% tidyr::gather(geneset, value, 36:ncol(ssgsea_plotdf01)) %>% group_by(ROIType, Timepoint, geneset) %>% summarise(AvgEnr = mean(value)) %>% tidyr::spread(geneset, AvgEnr) %>% as.data.frame()## `summarise()` has grouped output by 'ROIType', 'Timepoint'. You can override
## using the `.groups` argument.ssgsea_sdev = ssgsea_plotdf01 %>% tidyr::gather(geneset, value, 36:ncol(ssgsea_plotdf01)) %>% group_by(ROIType, Timepoint, geneset) %>% summarise(Sdev = sd(value)) %>% tidyr::spread(geneset, Sdev) %>% as.data.frame()## `summarise()` has grouped output by 'ROIType', 'Timepoint'. You can override
## using the `.groups` argument.rownames(ssgsea_avg) = paste(ssgsea_avg$Timepoint, ssgsea_avg$ROIType, sep="_")
avgexp_plot = ssgsea_avg[order(ssgsea_avg$ROIType, ssgsea_avg$Timepoint),c("ROIType", "Timepoint", names(geneset))]
rownames(avgexp_plot) = paste(avgexp_plot$Timepoint, avgexp_plot$ROIType, sep="_")
avgexp_plot0 = avgexp_plot[,3:ncol(avgexp_plot)]; rownames(avgexp_plot0) = rownames(avgexp_plot)
pheatmap::pheatmap(avgexp_plot0, #%>% scale(),
color = colorRampPalette(rev(brewer.pal(n = 10, name = "BrBG")))(100),
border_color = "grey60",
cluster_cols = T, cluster_rows = F,
legend = TRUE, scale = "column",
#annotation_row = annoCol_df %>% t(),
#annotation_colors = annoCol2,
show_colnames = T, show_rownames = T)
pheatmap::pheatmap(avgexp_plot0[c(1,3,5,7),]/avgexp_plot0[c(2,4,6,8),],
color = colorRampPalette(rev(brewer.pal(n = 10, name = "BrBG")))(100),
border_color = "grey60",
cluster_cols = TRUE, cluster_rows = FALSE,
scale = "column",
show_colnames = TRUE, show_rownames = TRUE)
avgexp_df_gg <- avgexp_plot0 %>% tidyr::gather()
avgexp_df_gg$ROIType <- factor(rep(avgexp$ROIType, length(unique(avgexp_df_gg$key))), levels=c("OE", "IE", "OD", "VA"))
avgexp_df_gg$Timepoint <- factor(rep(avgexp$Timepoint, length(unique(avgexp_df_gg$key))), levels=c("Pre", "Post"))
sdexp_df_gg <- ssgsea_sdev[,3:ncol(ssgsea_sdev)] %>% tidyr::gather()
avgexp_df_gg$sdev <- sdexp_df_gg$value
ggplot(avgexp_df_gg, aes(x=Timepoint, y=value, fill=Timepoint)) +
geom_bar(stat="identity", alpha=0.8, position=position_dodge(0.9)) +
facet_grid(key~ROIType, scales ="free") +
scale_color_manual(values=c("#377eb8", "#e41a1c")) + scale_fill_manual(values=c("#377eb8", "#e41a1c")) + theme_bw() +
geom_errorbar(aes(ymin=value, ymax=value+sdev), width=0.2, position=position_dodge(0.9))